
function [p q] = implicit_theta(dt,beta)
global U  h  g  p  q  nt dx x h1 p_global L;

dt2= 1/dt;
I = eye(2*x);
% Periodic Boundary Conditions
sol(1:x,1) =p;
sol(x+1:2*x,1) =q;
%p_global(:,1)=p;

for n=2:nt+1;   
    rhs = (dt2*I - (1-beta)*L)*sol;
    sol=(dt2*I+beta*L)\rhs;    
    p = sol(1:x,1);
    q= sol(x+1:2*x,1);
%p_global(:,n)=p;
%     if rem(n,100)==0
%         
%         %refreshdata(h1,'caller') % Evaluate p in the function workspace
%      %   drawnow
%     end
    
end
%display('Completed Successfully');